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ABSTRACT 

Context, About 2/3 of the Be stars present the so called V/R variations, a phenomenon characterized by the quasi-cyclic variation 
of the ratio between the violet and red emission peaks of the H I emission lines. These variations are generally explained by global 
oscillations in the circumstellar disc forming a one-armed spiral density pattern that precesses around the star with a period of a few 
years. 

Aims. In this paper we model, in a self-consistent way, polarimetric, photometric, spectrophotometric and interferometric observations 
of the classical Be star f Tauri. Our primary goal is to conduct a critical quantitative test of the global oscillation scenario. 
Methods. We have carried out detailed three-dimensional, NLTE radiative transfer calculations using the radiative transfer code 
HDUST. For the input for the code we have used the most up-to-date research on Be stars to include a physically realistic description 
for the central star and the circumstellar disc. We adopt a rotationally deformed, gravity darkened central star, surrounded by a disc 
whose unperturbed state is given by a steady-state viscous decretion disc model. We further assume that disc is in vertical hydrostatic 
equilibrium. 

Results. By adopting a viscous decretion disc model for f Tauri and a rigorous solution of the radiative transfer, we have obtained 
a very good fit of the time-average properties of the disc. This provides strong theoretical evidence that the viscous decretion disc 
model is the mechanism responsible for disc formation. With the global oscillation model we have successfully fitted spatially resolved 
VLTI/ AMBER observations and the temporal V/R variations of the Ho- and Bry lines. This result convincingly demonstrates that the 
oscillation pattern in the disc is a one-armed spiral. Possible model shortcomings, as well as suggestions for future improvements, are 
also discussed. 

Key words. Methods: numerical - Radiative Transfer - Stars: emission-line, Be - Stars: individual: ( Tau - Polarization - Techniques: 
interferometric 



1. Introduction 

( Tauri is a northern, nearby Be star that has drawn the attention 
of generations of observers and theoreticians. Its position on the 
sky makes it a convenient target for both northern and southern 
telescopes, and, as a result, this star has a very rich observational 
history. £ Tau has some distinctive features that make it an im- 
portant laboratory for testing theoretical ideas about processes 
associated with the Be phenomenon. First, it has a circumstellar 
(CS) disc that has shown little or no secular evolution in the past 
18 ye ars or so dStefl et ai1l2009l hereafter Paper I; Rivinius ~et al.l 
2006). Second, it shows a very stable Vy/Q cycle, with a quasi- 
period of about 1430 days (Paper I; IStefl et all [20071) . In Paper 

Send offprint requests to: A. C. Carciofi, e-mail: carciofi@usp.br 

1 The V/R term refers to the quasi-periodic variation of the ratio be- 
tween the violet and red emission peaks of circumstellar emission lines. 



I, a compilation of spectroscopic, spectrophotometric, polari- 
metric and spectropolarimetric observations of £ Tau was pre- 
sented, covering 3 complete V/R cycles from 1992 to the present. 
The data was complemented by the first spectro-interferometric 
observations of this star, made in December, 2006 with AMBER 
at VLTI. The reader is referred to Paper I for details about the 
observations and initial data analysis. 

In this paper we report a detailed modeling of the avail- 
able data. We adopt a physical model for the central star and 
the CS disc that is consistent with the most up-to-date research 
on Be stars, as described below. For the central star we adopt a 
model that includes rotational effects such a s geometrical de- 
formation (Domici ano de Souza et al.L 120031) and the latitudi- 



It is observed in £ Tau and numerous other Be stars. See Paper I and 
references therein for more details. 
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nal dependency of the stellar radiation dTownsend et al.L 12004). 
Rotational effects are very important in determining the local ra- 
diation field at a given point in the stellar environment, since 
they redirect part of the flux toward the polar regions. For the 
unperturbed state of th e CS disc we adopt a steady-state, viscous 
decretion disc mod el dLee et all Il99lt iPorterl 1 19991; lOkazakil 
l200lt ICarciofilstaTl I2006TT In this model, the material is ejected 
from the stellar equatorial surface, drifts outwards owing to vis- 
cous effects, and forms a geometrical ly thin disc with nearly 
Keplerian rotation (Rivi nius et aUl2006l) . 

In order to model both V/R variations of the spectral lines 
and the VLT1/ AMB ER observatio ns, we use the global disc os - 
cillation model of OkazakH d 19971) and lPapaloizou et ail (Il992l) . 
In this model, a one-armed m = 1 oscillation mode is superposed 
on the unperturbed state of the disc. The model is based on the 
fact that iTi—l oscillation modes are the only possible globa l 
modes in nearly Keplerian discs s uch as Be discs jKato, 1983). 
It was fi rst applied to Be stars b y lOkazakil dl99ll) and then re- 
vised by Papaloizou et alj d 1992b . 

Although the observed characteristics of the V/R varia- 
tions were explained qualitatively well by the global oscillation 
model, it was not fully satisfactory for several reasons. First, 
the oscillation period is sensitive to the subtle details of dis c 
structure as well as stellar parameters (Fift & Harmaned 12006). 
Given that there are usually several stellar and disc parameters 
that are only loosely constrained, this high sensitivity means that 
it is always possible to obtain the same period by different com- 
binations of parameters. For example, as far as the period is 
concerned, the same effect is produced by decreasing the disc 
temperature or stellar radius, or by increasing stellar mass or ro- 
tation velocity, thus making model predictions ambiguous. Since 
we use a variety of observations to constrain our models in this 
investigation, we are able to significantly reduce these ambigui- 
ties. Second, although the model predicts oscillation modes that 
precess in the direction of disc rotation (prograde modes), as 
suggested by observations (e.g. jTelting et al.L 1 19941) . the modes 
are less confined to the inner part of the disc than expected from 
the observations. This can be a problem for isolated Be stars, 
for which the disc extends to a large distance, but we expect that 
this is not so serious for truncated small discs such as those in bi- 
nary Be stars like £ Tau. As a matter of fact, we will demonstrate 
in later sections that a less confined mode better reproduces the 
observed data. 

In order to critically test whether the global disc oscillation 
model and the viscous Keplerian scenario are applicable for the 
CS disc of £ Tau, it is necessary to rigorously solve the cou- 
pled problems of the radiative transfer, radiative equilibrium and 
statistical equilibrium in non-local thermodynamic equilibrium 
(NLTE) conditions for the complex 3D geometries predicted by 
the global oscillation scenario. A 3D approach for the radiative 
transfer is important to account for geometrical effects such as 
escape of radiation towards the optically thin polar regions, or 
the shielding of part of the disc by a local density enhancement. 
W e use, for this purpose , the comput er co de HDUST de s cribed 
in ICarciofi & Biorkmanl d2008l 120061) and ICarciofi etail d2008l 
120041) 

In Sect.|2]we explain our main model assumptions. In Sect. [3] 
we outline the adopted computational procedure and discuss in 
details the fitting of the observations reported in Paper I. Finally, 
in Sect. [4] and E] we discuss our results and present a summary 
of our conclusions. 



2. Model Description 

The computer code HDUST is a fully 3D, non-local thermody- 
namic equilibrium (NLTE) code designed to solve the coupled 
problems of radiative transfer, radiative equilibrium and statis- 
tical equilibrium for arbitrary gas density and velocity distribu- 
tions. The NLTE Monte Carlo simulation performs a full spec- 
tral synthesis by emitting photons from a rotationally deformed 
and gravity darkened star. The star is divided in a number of lat- 
itude bins (typically 100), each with its effective temperature, 
gravity and a spect ral shape give n by the appropriate Kurucz 
model atmosphere (Kurucz, 1994). After emission by the star, 
each photon is followed as it travels through the envelope (where 
it may be scattered, or absorbed and reemitted, many times) un- 
til it escapes. During a simulation, whenever a photon scatters, 
it changes direction, Doppler shifts, and becomes partially po- 
larized. Similarly, whenever a photon is absorbed, it is not de- 
stroyed; it is reemitted locally with a new frequency and direc- 
tion determined by the local emissivity of the gas. HDUST in- 
cludes both continuum processes and spectral lines in the opac- 
ity and emissivity of the gas. Since photons are never destroyed 
(absorption is always followed by reemission of an equal energy 
photon packet), this procedure automatically enforces radiative 
equilibrium and conserves flux exactly. 

The interaction (absorption) of the photons with the gas 
provides a direct sampling of all the radiative rates, as well 
as the heating of the free electrons. Consequently, an iterative 
scheme is adopted in which the rate equations are solved at the 
end of each iteration to update the level populations and elec- 
tron temperature. The process proceeds until convergence of 
all state variables is rea c hed. The interested reader is referred 
to ICarciofi & Biorkmanl d2006l) for details of the Monte Carlo 
NLTE solution. 



2.1. Geometry 

The geometry of the problem is defined in Fig. 1 . The star rotates 
counterclockwise with the rotation axis parallel to the z direc- 
tion, and is located at the origin of the cartesian system (x,y,z). 
The CS disc is assumed to lie in the equatorial plane (xy plane). 
Since we are investigating non-axisymmetric precessing discs, 
we define the x-axis so that it precesses together with the disc. 
Therefore, both the x and y directions are time-dependent. 

We define the (x',y',z') system so that the observer is located 
at x' — oo and the plane of the sky is parallel to the y'z' plane. 
This system is obtained by rotating the xyz system first by (f> 
degrees around the z axis and then by 90 - i degrees around the 
y' axis. The angle i (0 < i < 180°) is the viewing angle and the 
angle <p (0 < <p < 360°) describes the time-dependent position of 
the x-axis. 

To complete the geometrical description of the system we 
specify the angle y that the projection of the rotation axis on the 
plane of the sky, z', makes with respect to celestial north. We 
adopt the usual convention that y is measured east from north 
and that y lies in the interval -180 < y < 180°. The angles y, <p 
and i are all free parameters. 

2.2. The Central Star 

Th e adopted stellar param eters are listed in TableQ] As discussed 
bv lRivinius et al.l d2006), determining the properties of the cen- 
tral star is diffi cult because of the large Vsin/ (320 km s _1 , 
lYangetallfT990h and the presence of an optically thick CS disc. 
The adopted parameters for the central star, together with the 
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x' (observer) 




Table 1. Model Parameters for £ Tau 



Fig. 1. Geometry of the problem. The star is located at the center 
of the xyz system, with the rotation axis aligned with the z axis. 
The disc is in the xy plane. The observer lies along the x' direc- 
tion, which is defined by the viewing angle i and the azimuthal 
angle <p, which is the angle between the x axis and the projection 
of x' onto the xy plane. The plane of the sky corresponds to the 
y'z' plane, and z' makes an angle y with respect to celestial north 



disc parameters (Sect. 12.3b , form a consistent set of parameters 
in the sense that, as we shall see below, they can reproduce well 
all the observational properties. However, it must be emphasized 
that many stellar parameters are somewhat uncertain. For in- 
stance, equally good fits are obtained by varying the physical 
size of the star or the photospheric temperatures by about 10%. 

Given the large value of Vsin/, we adopt a rotationally 
deformed and gravity darkened central star instead of a sim- 
ple spherical star. Assuming uniform rotation, the stellar rota- 
tion rate is given by a single parameter V/V cr it, i.e., the ratio 
between the equatorial velocity and the critical speed V cr i t = 
[2GM/(3Rp)] l/2 , where R p is the polar radius. We assume 
the stellar shape to be a spheroid, which is a reasonable ap- 
proxi mation for t h e sha pe of a rigidly rotating, sub-critical 
star (Fremat et al., 2005). Given the stellar rotation rate and 
the polar radius, the equatorial radius, R e , is determined from 
the Roche approxim ation for the stellar surface equipotentials 
dFremat et all 12005b . For the grav ity darkening we u se the stan- 
dard von Zeipel flux distribution (von Zeipel ]l924t) . according 
to which F(6) oc g e g(ff) oc T* s (ff), where g & g(0) and T e ff(9) are 
the effective gravity and temperature at stellar latitude 0. The 
stellar rotation rate was assumed to be V/Vcrit = 0.8. For the po- 
lar radius we have adopted the value R p = 5.9 Rq, typical for 
an ev olved main sequence star of spectral type Bl IV dHarmaned 
[19881: Paper I). 

The value of y cr ; t in Table Q] was determined from the data 
fitting in order to best reproduce the emission line profiles 
(Sect. |3.2| i. The value of the stellar mass was then calculated 
from y ci i t and R p . Both M and V cr jt we obtain agree with the 
values tabulated by Har manecl d 1988b . 



Parameter 


Value 


Type 


Ref." 


Stellar Parameters 


Spectral Type 


BUVe shell 


fixed 


2,3 


*p 


5.9 Rq 


fixed 


2,4 


R e 


1.1 Rq 


fixed^ 




T? 


25 000 K 


fixed 




T. 


17 950 K 


fixecF 




W Vcrit 


0.8 


fixed 




Vast 


530 km s _I 


free 




M 


11.3 M Q 


fixed^ 




L 


7400 Lq 


fixed 




T e ff 


19 370 K 


fixed f 




Disc Parameters 


So 


2.1 g cm 1 


free 




Po 


5.9 x 10" 11 gem" 3 


fixecF 




fl, 


130% 


fixed 




H 


0.208 Rq 


fixed* 




V/R Parameters 


P 


1429 d 


Fixed 


2 


T 


50414^ 


Fixed 


2 


fe 


0.006 ~ 


Free 




Weak line force 


5.74 x lO^nrARe) ' 1 


Free 




M 


lO^Moyr 1 


Free 




a 


0.4 


Free 




S 


0.95 


Free 




Geometrical Parameters 


i 


95° 


free 




y 


32° 


free 




<f>o 


280° 


free 




distance 


126 pc 


fixed' 


6 


Binary Parameters 


Porb 


132.97 d 


fixed 


7 


e 





fixed 


7 


a 


257% 


fixed 


7 


q 


0.121 


fixed 


7 



1) this 



" References: 
i Harmaned dl98 
Harmanec (198 

* R s was calculated from R 



work: 
5) ICranme 



2) Paper I; 31 lYang et 



i 2j_rape 
3 1 120050 : 



et~aU ( | 1990 () : 4) 
et all 



the Roche approximation. 
c The equatorial temperature, 



; 6) IPerrvman et all dl997h : 7) 
p and V/ Vcrit for a rigidly rotating star in 



T t , was calculated from V/V clit and T p 



for a rigidly rotating star in the Roche approximation, 
was calculated from V cr jt . 

e 7" e ff corresponds to the average effective temperature of a rotation- 
ally deformed, gravity darkened star. 

f po is given by E (27r)- 1/2 tf 1 (Eq.QOj. 

g Ho was calculated from Eq. [8] assuming a mean molecu- 
lar weight of 0.6 and an electron temperature of 18 000K after 
ICarciofi & Biorkmanl J2006li . 

* For To we have chosen the peak value of cycle I (see Paper I). 

' d was considered a free parameter within the 113 — 148 pc 
range, corresponding to the accuracy in the distance determination the 
Hipparcos satellite. 



2.3. Disc Structure 

We assume a steady-state viscous decretion disc, for which the 
surface density is given by (e.g jBiorkman & Carciofil 120051) 



«->-5s(3f)'W-'] • 



where 



<zr = (x z + y z ) 



2x1/2 



(1) 



(2) 
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is the dist ance from the center of the star, a is the viscosity pa- 
rameter of Sha kura & Sunvaevl(ll973l) . M is the stellar decretion 
rate and Rq is an arbitrary integration constant. The isothermal 
sound speed is c s = (khT) l ^ 2 (p,muy 1 ^ 2 , where kt, is the Boltzman 
constant, T is the gas kinetic temperature, yU is the mean molec- 
ular weight, and ntn is the mass of the hydrogen atom. 

The integration constant Rq is a free parameter related to the 
physical size of the dis c. For time-dependent models, such as 
those of Okazaki (2007), Rq grows with time and thus Rq is re- 
lated with the age of the disc. 

£ Tau is a well-known single line binary ( Harm aneH 1 19841) 
and it is usually assumed that the tidal forces from the secondary 
physically truncate the disc at a radius R t , correspond ing to the 
tidal radius of the system dWhitehurst & Kingl 1199 lb . We as- 
sume that the disc is sufficiently old that R t <K Rq, in which case 
Eq.[T]can be written in a simple parameterized form 



E(^7) = X R 2 e m- V2 [(Ro/v) 1 ' 2 - l] « So (^-) 



(3) 



where So is a constant that measures the density scale at the base 
of the disc, and is written as 



E = 



(4) 



In Eq. Q] the disc density scale is controlled both by Rq and 
M. For a truncated disc system Rq is unknown, because the in- 
formation about the disc age is destroyed by the tidal disruption 
of the outer disc by the secondary. Therefore, in such a system M 
cannot be determined observationally, unless some other infor- 
mation about the outflow is available. From the continuity equa- 
tion 



M = 2nml.(m)v c 



(5) 



we see that the mass decretion rate is related both to the surface 
density and the radial velocity, v w . Thus, in order to determine 
the decretion rate the observations must provide a measure of the 
outflow velocity. 

This latter quantity is difficult to obtain directly from the 
observations, because the outflow velocities in the inner disc 
are proba bly several orders of magnitude lower than the orbita l 
speeds dHanuschikl[l994l2000trWaters & Marlborough! [19941) . 
However, the outflow velocity can, in principle, be determined 
indirectly. As discussed below, the eigenvalues of the global os- 
cillation mode depend on the viscosity parameter a. This param- 
eter, in turn, sets the outflow velocity (this can be seen by substi- 
tuting Eq.[3]in Eq. [5J. Therefore, if we can sufficiently constrain 
the global oscillation eigenvalues (see Sect. 0), the value of a 
and, as a consequence, the decretion rate could be obtained. 

From the binary parameters of Table [1] we find that the 
Roche radius of the system is ^R OC he = 144 Rq. Her e we have 
used t he approximate formula for the Roche radius by Eg gletonl 
d!983l) . which is given by 



^Roche - A 



0.49^ 



-2/3 



0.69<r 2 / 3 +ln(l +q- l ^Y 



(6) 



where a is the semi-major axis and q — M2/M1 is the 
binary mass ratio. Since the tidal radius is approximately 
0.9j^R OC he in circular bina ries with small mass ratio (q £ 0.3, 
IWhitehurst & KingMl991l) . we assume that the disc around C, Tau 
is truncated at R t = 13O/?0. 



We also assume that the disc is in vertical hydrostatic equilib- 
rium. In this case , it can be shown that the den sity of isothermal 
discs is given by dBiorkman & Carciofil 120051) 



(7) 



p(m,z) = p (w)expl- — 

where po(vj) is the disc density at the midplane (z 
disc scale height is given by 

H{m) = H a [ — 



0), and the 



(8) 



where H = c s Vj t R e 

We can find the radial dependence of po by considering that 
S is the integral of p in the z direction 



(9) 



(10) 



/CO 
p(vj, z) dz = V27r//po(w) ■ 
CO 

From Eqs.[3]to|9]we have 
p{m,z)= — exp -^-r . 

Before describing the global oscillation model, a few words 
about the HDUST solution of the gas state variables are war- 
ranted. From Eqs.[3]and|7]we see that the solutions for both the 
vertical hydrostatic equilibrium and the radi al diffusion depends 
on the disc temperature. In a recent work. ICarciofi & Bjorkman 
(2008) studied the hydrodynamics of non-isothermal viscous 
Keplerian discs and found that the density structure of such 
discs can strongly deviate from the simple isothermal solu- 
tion of Eqs. Q] and [7] The dissimilarities between the two so- 
lutions are most marked for the inner disc (m ^ 10 R e ), from 
whence the polarization continuum and most of the line flux 
comes. Indeed, a comparison between the emergent spectrum 
for the self-consistent solution with isothermal models revealed 
important differences betw een the two (see, e.g., Fig. 6 of 
ICarciofi & BjorkmanLl2008l) . 

The models presented in this work have an isothermal den- 
sity structure but are not isothermal, because for the assumed 
density distribution HDUST calculates the full NLTE and radia- 
tive equilibrium problem, thereby providing the gas temperature 
as a function of position in the disc. ICarciofi & Bjorkman ( 2008) 
demonstrated that such mixed models are a much better approach 
to the problem than a purely isothermal model. The reason why 
an isothermal density structure was assumed in this work lies in 
the fact that at the moment only an isothermal solution for the 
global disc oscillations is available. Lifting this inconsistency 
between the calculated temperature distribution and the assumed 
density will be left for future work. 

To model the VI R variations and the VLTI/ AMBER observa- 
tions reported in Paper I, Eq.[3]must be mod ifie d according to the 
global oscillation model of OkazakU d 19971) andlPapaloizou ^taTI 
(Il992l) . In calculating the gravitational potential, we take into ac- 
count the quadrupole contribution due t o the rotational deforma - 
tion of the rapidly rotating Be primary (Papaloiz ou et al.Lll992l) . 
We a lso take into acc ount the azimuthally averaged tidal poten- 
tial dHirose & Osakil[l993l) . The potential is then given by 



GM 



1 +ko 



^crit 



2 / s-2 

m 1 



R e j +q a 



1 + 



, (ID 



where Q* is the angular rotation speed of the star, Q cl i t = 
2(3/?p)~' V cl it, and &2 the apsidal motion constant. In this equa- 
tion the first term is the point-mass potential of the Be star, the 
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second term is the quadrupole contribution, and the third term is 
the azimuthally averaged tidal potential. 

The radial distribution of the rotational angular velocity 
D,(m) is derived from the equation of motion in the radial di- 
rection, and is written explicitly as 

+ ^(«f_,(°)*l" 2 , (12) 

d\n m \m) \R e 



under the approximation z 2 /m 2 <k 1, where we have included a 
hypothetical radiative force in the form of 



GM lm 

frad = — hr 



(13) 



where r\ and e are parameters cha racterizing the force due to a n 
ensemble of optically thin lines dChen & M arlborough, 1994). 
Then the associated local epicyclic frequency k(tu) is written as 



2Q 2Q + B7- 



dQ. 

T 

dm 



1/2 



1 fe (fi cnt ) Ue) 2q (a) 
(H\ 2 lm 



dlnE d 2 \nl. 
2— + 



d In m dlnm 



1/2 



.(14) 



The perturbed surface density is obtained by superposing 
a linear m — 1 perturbation on the above unperturbed state 
(Eq. [3} in the form of normal mode of frequency co that varies 
as exp[/(wf - <f>)]. For simplicity, the perturbation is taken to be 
isothermal. The linearized perturbed equations are then obtained 
as follows, 



i{cj - Q) + v n 



dm 



£' 1 d iVj. 



m"E dm 



., „, dv m d 
i(co-Ll) + —— + v m — 
dm dm 



v' m -2Clv', = 0,(16) 



; d \ Jf 
m dm ] 2 



2D. 



v m d 
i(co-D) + — +v m — 
m dm 



0. 



(17) 



where U is the Eulerian surface-density perturbation, and 
(v' m , v'^) is the vertically averaged velocity field associated with 
the perturbation. We solve Eqs. (TT3T>-(fT7b with a rigid wall 
boundary condition, {v' m , v'J = 0, at the inner edge of the disc 
and a free boundary condition, Ap = 0, at the outer disc radius, 
where Ap is the Lagrangian perturbation of pressure. 

It is instructive to consider the effects of some parameters on 
the m = 1 oscillations before solving the perturbation equations. 
For simplicity, we first neglect viscous terms. Then, the local 
dispersion relation is obtained from Eqs. (fT3ll-(fT7]i as 



(co - D - kv m ) L - n 1 = elk 1 , 



(18) 



where k is the radial wave number. In this case, the propagation 
region, where A; is a real number, is given as the region where 



(CO — Q) 2 - K 1 + V^K 1 Icl > (e.g., |Okazaki, 2000), which is 
explicitly written as 



to < Q — K 



1-1 — 



1/2 



(GM\ 1/2 



1/rflnZ t/ 2 lnZW//\ 2 3 
2 [dlnm + dlnm 2 j\m) + 4 q 



2 \n c n t ) \m) 2 ' \R e ) 2\cJ 



(19) 



Here, the first, second, and third terms are due to the pressure 
gradient force, tidal force and the quadrupole contribution of the 
potential around the deformed star, respectively. The fourth term 
is due to the radiation force and the last term is the contribu- 
tion due to advection. Near the star, the quadrupole term and 
the radiation term play the major role. Thus, the eigenfrequency 
is mainly determined by these terms and the pressure gradient 
term. The larger the contribution of the quadrupole term and the 
radiation term and/or the smaller the sound speed (or tempera- 
ture), the higher the eigenfrequency. On the other hand, the other 
terms, i.e., the tidal term and the advective term, as well as the 
pressure gradient term affect the oscillation behavior in the outer 
part of the disc. The larger the contribution of these terms, the 
less the eigenmode is confined, and vice versa. 

Next, we consider the effects of viscosity. When the mass- 
decretion rate M and the surface density at the inner disc ra- 
dius Eq are fixed, the radial velocity v m is proportional to the 
viscosity parameter a. Thus, one effect of viscosity is to make 
the o scillation mode le ss confined, in terms of the advective 
term. iKato et all ([T978) showed that m = 1 modes in nearly 
Keplerian discs, which are neutral if the viscosity is neglected, 
become overstable when viscous effects are t aken into account . 
The g rowth rate is proportional to a (see also Negu eruela et all 
I2001I) . 

It is important to note that our model provides prograde 
modes for a plausible range of parameters. The fundamental 
mode, in general, has the eigenfrequency of the order Q(/? e ) - 
k(R & ). From Eq. (fT9l it is obvious that Q.(R e ) - K(R e ) is positive 
even if no radiation force is included. Thus, there is no way to re- 
produce the observed V/R period of £ Tau by retrograde modes, 
which have negative eigenfrequencies. 

It should also be noted that linear models like ours cannot 
predict the amplitude of the oscillations. Thus, for the purpose of 
comparison with various observations, we will assume that the 
nonlinear perturbation patterns are similar to the linear eigen- 
functions obtained from the above equations, and take the max- 
imum value of the perturbed part of the surface density, to 
be a free parameter, 6. 

3. Results 

3. 1 . Modeling Procedure 

The global oscillation model described above is not a trivial 
problem for a radiative transfer solver. A very large number of 
grid cells (typically about 100000) is needed to accurately de- 
scribe the m, z and dependence of the density, of the gas state 
variables and of the radiation field (for details about the adopted 
cell structure in HDUST, suc h as radial and latitude spacing, see 
ICarciofi & Bjorkmani |2006). Solving the radiative equilibrium 
and NLTE statistical equilibrium problems for each cell involves 
a long series of iterations to accurately determine the coupling 
between radiation and state variables in the dense environment 
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of the inner disc. This, and the subsequent calculation of the ob- 
servables, takes between 2 and 3 days for a single model in a 
cluster of 50 Pentium IV computers. 

This is prohibitively long for the iterative procedure neces- 
sary to fit the observations. For this reason, we have initially 
modeled the global disc properties using a simple 2D viscous de- 
cretion model using Eqs. (01 and (0. This 2D analysis allowed 
us to place initial constraints on several parameters of the sys- 
tem, including So, V C rit, h and d (Sect. [3T2b . 

Once a suitable 2D model was found, we used it as a start- 
ing point for the detailed 3D modeling described in Sect. 12.31 
With this 3D model we performed a simultaneous analysis of 
the VLTI/ AMBER interferometry and the V/R properties of the 
H I lines, that allowed us to constrain well the global disc oscil- 
lations parameters (Sect. 13.31 ). 

3.2. Two-Dimensional Viscous Decretion Disc Model 

As discussed in Paper I, £ Tau has shown little or no secular 
evolution since the onset of the current V/R activity started in 
1992, and it is argued that the global disc properties, basically 
its density scale, have been approximately constant in the period. 
Therefore, modeling the data with a 2D viscous decretion disc is 
per se a useful exercise that can provide insight into the average 
properties of this disc. 

Initially, a set of observations that are representative of the 
period from 1992 through 2008 must be defined. For both the 
spectral energy distribution (SED) and polarization in the 0.3 
— 1 /im region, we used the PBO spectropolari metric observa- 
tions made in March, 1996 (IWood et all ll997l) . The flux and 
polarization levels were scaled in order to match the average V- 
band flux and polarization of the period ((V) = 3.02 + 0.07, 
(p v ) = 1.46 ± 0.08%, see Paper I). To extend the SED to the 
IR, we used archival data from the IRAS point s ource cata- 
log dBeichman et all [19881) and the 2MASS catalog dCutri et all 
120031) . 

Because the line profiles are variable and strongly dependent 
on the disc asymmetries caused by the global oscillations, we 
have fitted only the average Ha peak separation in this period 
(240km s 1 , Paper I) and the average Ha EW (-15. 5 ±0.1 A, Paper 
I). 

One of the most beautiful characteristics of the viscous de- 
cretion disc model is its relatively small number of parameters. 
To describe the steady-state structure of the disc we need to spec- 
ify the decretion rate, M, the viscosity parameter, a, and the age 
of the disc that is enclosed in the parameter Rq (Eq.[T|i. We fur- 
ther need to specify the stellar critical velocity, V cr i t , which is a 
measure of the stellar mass and affects the disc structure by set- 
ting its scaleheight (Eq.[8]l and rotation speeds. In the case of the 
disc around ( Tau we assume, as explained in Sect. 12.31 that the 
disc is truncated by a companion star. This mechanism erases the 
information about the disc age and makes the problem undeter- 
mined. This basically means that M, a and Rq are all contained in 
a single parameter, So, which determines the disc density scale. 
Therefore the disc is described by only two free parameters, So 
and V cl i t . In addition to these two physical parameters we must 
specify two geometrical parameters, the viewing angle, i, and the 
distance to the star, d. 

For this analysis several tens of models were run, covering 
a large range of values for So, V C rit and i. For each model, the 
value of d was computed in order to match the observed flux lev- 
els. Since d is known to be in the range between 113 — 148 pc, 
which corresponds to the accuracy in the distance determination 
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Fig. 2. Emergent spectrum of ( Tau. The dark grey lines and 
symbols are the observations and the black lin es represent t he 2D 
model results. Top: Visible SED (data from I Wood et al.Lll997l 
scaled in flux to match the average V-band magnitude from 1992 
to the present). Middle: IR SED (data from 2MASS - squares, 
and IRAS - circles ). Bottom: Continuum polarization (data from 
IWood et all Il997l also scaled in level to match the average V- 
band polarization from 1992 to the present). In the two upper 
panels, the light grey line corresponds to the unattenuated stellar 
SED 



by the Hipparcos satellite (iPerrvman et all 1 1997b . models that 
did not produce d in this range were discarded. 

Figs. |2] shows our best fit to the data using the 2D viscous 
decretion model. This best fit was achieved for the following 
parameters 

S = 1.7 gem" 2 , 
i = 85° , 

y crit = 53010ns- 1 , 
d = 126 pc, 
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Fig. 3. Continuum polarization for £ Tau. The model results for 
three different viewing angles, as indicated, are compared with 
the observations 



in addition to the fixed stellar parameters listed in the first part 
of Table Q] We recall that for all models the disc was truncated 
at 130%. 

The above value for Zo corresponds to po = 5.9 x 
1CT 11 g cirT 3 . The number density of particles can be estimated 
assuming a mean molecular weight of 0.6, typical for ionized 
gas with solar chemical composition (recall that our models have 
only H). We obtain a value of 5.9 x 10 13 cirT 3 , whi ch is character- 
istic o f Be stars with dense discs, such as 6 Sco dCarciofi et al.L 
120061) . 

Before discussing our results, it is useful to recall the physi- 
cal processes that control the emergent spectrum. In the visible, 
the SED depends on the spectral shape of the stellar radiation, 
which, for a rotating star, is controlled by the latitude-dependent 
effective temperature. The stellar radiation is modified by several 
processes in the CS disc: electron scattering, bound-free absorp- 
tion and free-bound emission by H I atoms, free-free absorption 
and emission by free electrons, and line absorption and emission 
by H and other elements. The absorption and emission processes 
depend on the atomic occupation numbers and electron tempera- 
ture, which, in turn, depend on the radiation field in a non-linear 
and complex way. 

In the IR the SED is controlled mainly by the free-free emis- 
sion of the free electrons. The continuum polarization is pro- 
duced by electron scattering, but is also modified by the absorp- 
tion and emission processes in the disc. 

The absorption, emission and scattering of the radiation de- 
pends also on the geometry and velocity structure of the CS 
material. For instance, the linear polarization is strongly depen- 
dent on the number density of electrons, the disc flaring and 
scaleheight. The slope of the IR SED, on the other hand, de- 
pends mainly on the radial profile of the density, but also on the 
disc fl aring and temperature distribution dCarciofi & Biorkmani 
120061) . 

In view of the complex dependence of the emergent radia- 
tion on the stellar and disc parameters, the overall good match 
between the model and observations is remarkable. The model 
reproduces the SED longward of the Balmer jump (3646 A) all 
the way down to the longest IRAS wavelength. The only part of 
the SED which is not well reproduced is the size of the Balmer 
jump. In our models, the size of the jump is too large, meaning 
that either the amount of neutral hydrogen in the line of sight is 



a little overestimated or the geometry of inner disc is not well 
described in our models (see Sect. |4]i. 

The predicted polarization matches the observations well. 
The size of the Balmer jump and the level and slope of the 
polarization in the Paschen continuum show good agreement. 
However, in the Brackett continuum our models overestimate the 
polarization by approximately 10%. 

We obtain an inclination angle of i = 85° for our best-fit 
model. The fact that C T au is a well-known shell star (Porter, 
1996; Rivinius 'et all l2QQ6h supports our findings. In Fig. [3] we 
compare the model polarization for 3 different inclination an- 
gles, 83°, 85°and 87°. Differences in i as small as 2° produce 
significant changes in the polarization, which indicates that the 
inclination angle is well constrained by our analysis. 

The extent that the CS disc modifies the stellar radiation can 
be assessed by comparing the model SED with the unattenu- 
ated stellar radiation in Fig. [2] Since the star is viewed close to 
edge on, this causes extinction of UV and visible radiation that 
is reemitted at longer wavelengths. Most of the reprocessed ra- 
diation escapes in the polar direction, due to the fact that the ver- 
tical optical depth is much smaller than the radial optical depth. 
However, a portion of this radiation will escape radially and will 
contribute to the large IR excess observed 3 mag at 25 fim). 

The Ha peak separation and EW for the best-fit model are 
240km s _1 and -14 A, respectively, in good agreement with the 
average values for this time period. We have chosen V cr j t as a 
free parameter and, therefore, no a priori assumption about the 
stellar mass have been made. V cr j t is well constrained by the 
model, as both the CS emission lines and the continuum po- 
larization depend on this value (recall that V cr jt affects the disc 
scaleheight). V cr i t and the assumed R e - 7.7 Rp, giv e a stellar 
mass of 1 1.3 Mq, in agreement with Harmanec ( 1984). This fur- 
ther supports the values of the parameters adopted for ( Tau. 

The good fit that was obtained using a relatively simple 
model provides strong observational support to our assumptions. 
The power-law description for the surface density proved suc- 
cessful in determining the slope of the IR SED, and the hydro- 
static equilibrium assumption successfully determined the cor- 
rect shape and level of the continuum polarization. Therefore, 
we conclude that a steady-state viscous decretion disc, truncated 
at R t = 130 Rq, is a good description for the average properties 
of the CS disc around £ Tau. 

3.3. Global Oscillation Model 

We now present our analysis of the VLTI/ AMBER observations 
and the V/R properties of £ Tau using the global oscillation 
model. In addition to the parameters Xo, Vcrit, d, and i, that we 
have discussed above, we now consider the parameters kz, 5, M, 
a and the weak line force, so that we can consider the m = 1 
density perturbation within the disc. Also, the time-dependent 
position of the m = 1 perturbation pattern with respect to the 
direction of the observer, described by the parameter <p, must be 
determined (see Fig.Q]). 

Different aspects of our solution are illustrated in Figs. [4] to 
[9] and the parameters of the best fit obtained are summarized 
in Table Q] Before continuing to describe our results, some ex- 
planations about our modeling procedure are warranted. For a 
given model, HDUST computes the emergent spectrum for an 
observer whose position is specified by i and (f>. To simulate the 
temporal dependence of the observables as the density wave pre- 
cesses around the star, we compute, for a given i, the emergent 
spectrum for several values of <p. The global oscillation parame- 
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Fig. 4. log(V/R) vs. time for the Ha (black), Bry (green) and 
Brl5 (red) lines. The points correspond to the observations, and 
the lines to our best fitting model. The red dashed line corre- 
sponds to the Brl5 V/R curve artificially shifted in phase by -0.2 
to match the observations 

ters used in the simulation are then evaluated by comparing the 
theoretical amplitude and shape of the V/R cycle to the observa- 
tions. 

The time-dependent position of the density wave is related to 
the V/R phase as follows. According to the convention of Paper I, 
we define the V/R phase t such that t = at the V/R maximum. 
The ephemerides of the V/R cycle, given by the period, P, and 
the modified Julian date for r = 0, To, were determined in Paper 
I and are given in TableQ] We introduce the parameter 0o, which 
gives the position of an arbitrary point of the spiral pattern at 
t = 0. We chose this point to be the minimum of the density 
perturbation pattern at the base of the disc. The angle </> is related 
to the phase by 



(Note that the angle (f> decreases with time for the prograde os- 
cillation mode assumed here). Since the phase is related to the 
time by r = (T - T$)P~ l , we find the following expression for 
the temporal dependence of (/> 

^ = 27r|l-I_I^ + 0o. (21) 

This equation gives the position of the minimum of the density 
pattern at the base of the star as a function of time. For a given 
model, a value of <po is obtained by fitting the V/R phase and the 
interferometric observations. 

3.3.1. Constraints from the V/R Observations 

In Fig. 2] we show our best fit to the V/R cycle of the Ha, Bry 
and Brl5 lines. For Ha, the model reproduces the first part of 
the V/R cycle (0 < t < 0.5), including the shape of the curve 
and the observed V/R maxima and minima. In the second half of 
the cycle, between 0.5 < t < 0.8, the observed V/R ratio departs 
from the quasi-sinusoidal shape that is characteristic of the first 
half (Paper I), and this is not reproduced by the model. 

The behavior of the V/R properties of the IR lines are similar 
to that of Ha, but the smaller number of observations make the 
results less conclusive. Unfortunately, there are only three ob- 
servations in the first half of the cycle, but their V/R ratios agree 



with the predicted values (see points for MJD w 53 700). The 
second half of the cycle is better covered. For Bry, the model 
seems to reproduce the data, but the V/R values tend to cluster 
above the predicted curve, with the exception of the two points 
for r > 0.8. For Brl5, all the data points are above the predicted 
values in the second half of the cycle. 

One interesting feature of the observations that is clearly 
seen in Fig. [4] are the phase differences between the cycles of 
different lines. Th ese phase differences ha ve been already ob- 
served in £ Tau by Wis niewski et alj d2007f) . who suggested that 
they might be associated with the helicity of the global oscilla- 
tion. The idea behind this suggestion is that the more optically 
thin infrared lines form closer to the star in the disc, whereas 
the Ha line forming region extends further out into the disc. 
Because of the strong helicity of the density perturbation pattern 
(see, e.g., Fig. |5]), the average azimuthal morphology of the in- 
ner disc is quite different than the morphology of the outer disc. 
Therefore, the observed phase differences would be a result of 
the different average morphology of each line formation site. 

The phase differences between the Ha and Bry lines are well 
reproduced by our models, which provides convincing evidence 
in support to a spiral morphology of the density perturbation pat- 
tern in the disc of £ Tau. The phase of the Brl5 line is not well 
reproduced, however. This line is formed in the very inner part 
of the disc and it seems, therefore, that the global oscillation 
model does not describe correctly the shape of the spiral pattern 
in this region. In Fig. [4] the dashed red line shows the model V/R 
curve of the Brl5 line shifted in phase by -0.2. This artificially 
displaced curve reproduces the observations. This may indicate 
that the model predicts the correct amplitude of the oscillations 
in the inner disc, even though the shape, which sets the phase, is 
not correct. 

As mentioned above, the global oscillation model failed to 
reproduce the Ha V/R cycle for 0.5 < t < 0.8. However, this 
may not be a problem with the global oscillation model, since 
the 0.5 < r < 0.8 phase interval is where the complex triple- 
peaked Ha profiles become more prominent. In Paper I, we have 
presented evidence that the triple-peak anomaly arises from the 
outermost part of the disc, and that it is probably a line-of-sight 
effect owing to the large inclination angle of the system. The 
failure to reproduce this part of the V/R cycle is likely due to the 
fact that a physical process is missing in our analysis. 

Fig. [5] illustrates the properties of our solution for the disc 
global oscillation. The oscillation mode is subject to two very 
stringent constraints: it must have a precessing period equal to 
the observed V/R period (1 429 days, see Paper I) and it must be 
able to reproduce the observed amplitude of the V/R variations. 
This latter requirement is associated with the degree of confine- 
ment of the mode, which essentially means how far out into the 
disc the perturbation goes. In Fig.|5]we compare the theoretical 
V/R curve for two modes with the same precession period but 
different degrees of confinement. The model shown in the top 
left is the best fit model of Table Q] This model corresponds to 
a mode with a small degree of confinement, in which the os- 
cillations extend all the way down to the outer rim of the disc. 
The other model is shown as an example of a mode that also has 
the same precession period but it is much more confined. From 
the comparison of the predicted V/R ratios for each mode (bot- 
tom panel of Fig. [5j we see that this latter mode can be readily 
discarded on the basis that it cannot account for the observed 
amplitude of the V/R variations. 
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Fig. 5. Comparison of the results for modes with different de- 
grees of confinement. The top left plot shows the density per- 
turbation pattern for the V/R parameters of Table [1] The pic- 
ture depicts the disc as seen from above and direction of preces- 
sion of the density wave is counter-clockwise. A more confined 
mode, with a weak line force of 0.04(-nr//? e ps) > a = 0.3 and 
ki = 0.006, is depicted in the top right panel. The bottom panel 
shows the predicted amplitude of V/R ratio for each model (solid 
line: less confined model; dotted line: more confined model). 



3.3.2. Contraints from the AMBER/VLTI Observations 

Perhaps the strongest evidence in support of the global oscilla- 
tion scenario comes from the very good fit of the interferomet- 
ric data. The formal fit of the observed visibilities and phases is 
shown in Fig.[6]for three baselines (marked as A, B and C, see 
Fig. [7] for their orientation with respect to the disc). In the fig- 
ure, the left column shows the phases and the middle column the 
square visibilities for each baseline. Because of the poor abso- 
lute calibration of our dataset, both observables were normalized 
to the value obtained in the continuum. The difference in am- 
plitude of the different profiles are explained by both the base- 
line lengths (93m, 53m and 130m, respectively) and orientation 
(43°, 99° and 63° East from North, respectively). For instance, 
baseline A corresponds to an orientation close to the minimum 
elongation axis of the disc, which explains the small drop in vis- 
ibility. The quality of the fit is quantitatively expressed by the x 2 
map shown in the right column. 

For illustration purposes we have converted the three in- 
terferometric phases into an astrometric displacement for each 
spectral bin across the line (top-right panel of Fig. |6j. Open and 
filled symbols indicate the displacement along the disc at min- 
imum and maximum elongation axes, respectively. The typical 
S-shape astrometric signature of a rotating disc is clearly dis- 
torted: the redshifted part of the disc extends further away from 
the central star than the blueshifted part (300 //as and 150 yt/as re- 
spectively). This is a direct, illustrative proof of the presence of 
an asymmetry in the disc. 



The fits of the visibilities and phases were carried out in the 
following way. For a given disc model and inclination angle, we 
computed synthetic images for several wavelengths around the 
Bry line (see Fig. [7]) for several values of <p, corresponding to dif- 
ferent V/R phases. This set of synthetic images were then rotated 
for several values of y, corresponding to different orientations of 
the disc in the skyplane. Next, the visibilities and phases for each 
of the modified images were computed and compared to the ob- 
served visibilities and phases. A% 2 minimization procedure was 
carried out that produced the values of y and that best fitted 
the datfl Finally, we manually explored the parameter d in the 
range 1 13 - 148 pc. However, we have always found that this had 
only a marginal impact on our fit and, therefore, we have chosen 
to use the distance determined from the 2D analysis, d — 126 pc. 

This procedure was repeated for several disc models and in- 
clination angles. An analysis of the flux, polarization and V/R 
properties was also carried out for the same models. The model 
parameters of Table Q] correspond to the values that best repro- 
duced all data analyzed, including the interferometry. 

Our data analysis indicates that the interferometry is much 
more sensivite to some model parameters than others. The most 
sensitive parameter, and hence the best constrained, is undoubt- 
edly the disc orientation y — 32 + 5°. As discussed in Paper I, 
this value agrees well with previous polarimetric and interfero- 
metric determinations. In particular, it is compatible with the re- 
cently published va lue y = 37 + 2°, o btained in the K' band with 
the CHARA array (iGies et al.U2007l) . Additionally, as shown in 
the middle-right panel of Fig. [6] the overall size o f our models 
also agrees with that measured bv lGies et al.l(l2007l) . We note the 
CHARA and VL TI/AMBER estimates are independent because 
IGies et all d2007l) used absolutely calibrated broad band visibil- 
ities, whereas we have used differential phases and visibilities 
across a line. 

The x 1 ma P °f Fig. [6] shows that the parameter <p is not as 
well constrained. From the fit we find that the position of the 
density wave at the time of the observation is <p(T = 54 081.3) = 
90 + 25°. The large uncertainty comes from the combined ef- 
fects of the high inclination of the disc and the fact that the 
density perturbations are spread over a large range in azimuth 
(see Fig. [7]). Using Eq.|2T|we find that the position of the density 
wave at t = 0, as given by the interferometry, is (p™ 1 = 295 ±25°. 
Despite this large uncertainly, this value is in agreement with our 
best-fit value of 0o = 280° obtained primarily from the fit of the 
V/R phases. We emphasize that spectroscopic and interferomet- 
ric determinations of the 0o are completely independent, making 
the above agreement relevant. 

The analysis of the interferometric data for different values 
of i (not shown) largely confirms the value for the inclination 
angle obtained from the 2D model (Sect. 13.21) . In addition, in- 
terferometry has allowed us to lift a historical degeneracy in the 
determination of i. Recall from Fig. Q] that the inclination angle 
goes from to 180°. For i < 90° the top of the disc is facing the 
observer and the contrary is true for i > 90°. 

In the case of an axisymmetric disc, it is really difficult to 
distinguish between i > 90° and i < 90°. For an infinitely thin 
disc this distinction is impossible, since the disc aspect would 
be identical for the two cases. For a disc of non-zero geometri- 
cal thickness, radiative transfer effects create subtle differences 
in the aspect of each disc, but the differences are probably too 



2 For the x 2 calculation only the spectral channels in the range 
[-700, +700] kms~' were used. The value of^ 2 depends on the num- 
ber of spectral channels, and in the continuum these channels are always 
well matched and therefore tend to artificially reduce the^ 2 . 
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Fig. 6. Fitting of the interferometric observations for the spectral region around the Bry line. In all panels the solid line represents 
our best-fit model of Table [1] and the dashed line a model with i - 85°. Left and Middle: Differential phases and differential 
visibilities obtained with AMBER/VLTI, for three interferometric baselines. The baseline orientations are indicated in Fig. [7] Top- 
Right: Interferometric phases converted into astrometric shift (photocenter displacement) vs. wavelength. Open and filled symbols 
represent displacement along the disc minimum and maximum elongatio n axis respective ly. Middle-Right: Angular diameter derived 
from a Gaussian fit of the CHARA K' band continuum observations of iGies et al .1 (120071) . Bottom-Right: x 2 map showing the fit for 
the y and <p angles for the i — 95° model, with contours plotted for^ 2 = 4 and^- 2 = 6 



small to be detected. However, for a non-axisymmetric disc this 
distinction is possible, provided two conditions are fulfilled: 1) 
spatially and spectrally resolved observations are available in or- 
der to determine the rotation direction of the disc in the plane of 
the sky, so that the northward direction can be defined, and 2) 
the disc asymmetry is chiral, i.e., the disc is not superimposable 
on its mirror image. We note that we adopt the usual convention 
that the angular momentum vector points north. 

Both of the above conditions are fulfilled for f Tau. The 
VLTI/AMBER observations unambiguously determine that the 
western side of the disc approaches the Earth. Thus, according 
to the definition of north and the value of y derived from the in- 
terferometry, we know that the north of the system is oriented 
32° east of celestial north. Condition 2 is also fulfilled, since 
the disc is highly non-axisymmetric because of the spiral den- 
sity perturbation produced by the global m = 1 oscillation. From 
the detailed model fitting of the interferometric observations, we 
find that the model with viewing angle i = 95° (x 2 = 1 -4) fits the 
observations better than the corresponding model with i - 85° 
(X 2 = 2.4, see also significant differences between the models in 
Fig.|S). Therefore, the VLTI/AMBER data allow us to determine 



that the bottom part of the disc of ( Tau faces the Earth. Note that 
this degeneracy is nearly impossible to be lifted with broadband 
interferometry data only, such as the CHARA data shown in the 
right-middle panel of Fig. [6] We believe that this is the first time 
the degeneracy in the determination of ; was lifted for a Be star. 

3.3.3. Geometrical Considerations 

The geometrical properties of ( Tau derived above are illustrated 
in Fig. Q In the top left panel we plot the density perturbation 
pattern, as seen from above the disc. The figure marks the posi- 
tion of the observer's line of sight for 5 V/R phases. The t = 0.57 
line of sight corresponds to the time of the interferometric ob- 
servations. The top center panel is a sketch of the disc density 
projected onto the skyplane at the time of the VLTI/AMBER 
observations. The synthetic image for the IR continuum (top 
right panel) mimics the projected disc density as expected. The 
bottom panels of Fig. [7] show model images for different wave- 
lengths across the Bry line. The combination of the disc asym- 
metries and the effects of rotation generates a very complex 
brightness distribution. 
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Fig. 7. Top left: Density perturbation pattern of our best-fit model, as seen from above the disc. The line of sight towards the observer 
is marked for 5 V/R phases. The VLTI/ AMBER observations corresponds to r = 0.57. The two circles mark the position on the 
disc for m — 2R e and m — 5 R e . Top center: Density perturbation pattern projected onto the plane of the sky, at the time of the 
VLTI/ AMBER observations. The solid line indicates the position of the stellar rotation axis. The dashed lines marks the orientations 
of the 3 baselines shown in Fig. [6] Top right: Continuum synthetic image at 2.16 pm. Apparent from the plot is the brighter southern 
hemisphere of the star, which is little affected by the presence of the geometrically thin disc. Bottom panels: Synthetic images for 
three wavelengths across Bry, as indicated 



From Fig. [7] we see that in December 2006 (r = 0.57) the 
base of the spiral arm was behind the star, but the bulk of the 
material was located in the north-east. This is consistent with 
the V/R < 1 value at the period of VLTI/ AMBER observations 
(recall that the western side of the disc is approaching us). 



3.3.4. Polarization and Flux 

Besides the line V/R variations, the model also predicts flux and 
polarization variations across one V/R cycle. This can be readily 
seen from Fig. [8] where we compare the observed SED, polar- 
ization spectrum and Ha line profile with the model predictions 
for different phases. The polarization variations are further illus- 
trated in Fig. [9] 

What causes the variations in polarization is the occultation 
of part of the disc by the star in this nearly edge on system. Our 
models show that the polarization comes mainly from the region 
of the disc between m = 1 to uj « 5 R e , From the diagram in 
Fig.|7J we see that the bulk of the density between those radii is 
in front of the star for r w 0.5. Since polarization is roughly pro- 
portional to the number of scatterers, the polarization is expected 
to be greatest at this phase, as in indeed the case (Fig. [9). As the 
density perturbation pattern precesses, eventually the bulk of the 
material will move behind the star at t » 0. The light scattered 
by this material will be partially absorbed by the star causing a 
decrease in the polarization. 



In Fig.|9]the predicted modulation of the polarization is com- 
pared to the observations (data from Paper I). We have folded 
the data into phase using the V/R ephemerides of Table Q] We 
plot only observations from 1992 to 2004. The polarization data 
seems to present a big challenge to the model. The model pre- 
dicts large polarization variations of up to 40%, and this is not 
observed. We will come back to this important point in Sect. [4] 

The flux excess in visible wavelengths comes from light 
emitted by the recombination of a free electron with a proton, 
commonly called bound-free radiation. The recombination rates 
are roughly proportional to the square of the density and fall very 
rapidly with distance from the star. Thus the bulk of excess flux 
in the visible comes from a very small region of the disc, out to 
approximately 2 stellar radii. If we consider this small region, 
the density enhancement due to the global oscillations is behind 
the star for r » 0.5 and in front of it for r w 0, just the opposite 
situation described above for the polarization (Fig. [7]). This ex- 
plains why the flux is smaller for r = 0.58 in Fig. [8] The model 
does not predict important variations of the IR fluxes (second 
panel of Fig.[8]l because in this case the emission comes from a 
much larger fraction of the disc and occultation by the star has a 
smaller effect. 

Unfortunately, to our knowledge no systematic photometric 
monitoring was carried out for £ Tau after the onset of the V/R 
variations in 1992, so our predictions for the brightness varia- 
tions cannot be tested at the moment. 
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Fig. 8. Emergent spectrum from f Tau. The observations are 
shown as the dark grey lines and the V/R phase of each ob- 
servation, r bs, is indicated in each panel. Also plotted are the 
model predictions for three different V/R ph ases, as indicated i n 
the top panel. Top: Visible SED (data from I Wood eTafl I1997I) . 
Second from top: IR SED (data from the 2MASS and IRAS cat- 
alogs. Second from bottom: Continuum polarization (data from 
IWood et aUll997l) . Bottom: Ha line profile (data from Paper I). 
As in Fig. [2] the light grey line corresponds to the unattenuated 
stellar SED 
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Fig. 9. Temporal evolution of the polarization across the V/R 
phase. The model polarization (line) is compared with the ob- 
servations (points with error bars) 



We end this section with a discussion about the role of the 
parameter So in our 3D analysis. When the global oscillations are 
considered, So represents the azimuthally averaged density scale 
of the disc. We found from our 3D analysis a value of So 25% 
larger than that obtained from the 2D analysis (see Sect. |3.2b . 



4. Discussion 

Our analy sis of the data u sin g the global oscillation sc enario pro- 
posed by Okaz akld 1997b and lPapaloizou et all d 19921) shows that 
this model has both strong and weak points. We have been able 
to reproduce many observational characteristics of £ Tau. Of par- 
ticular relevance is the overall good fit of the V/R cycles of Ha 
and Bry lines and the excellent fit to the VLTI/ AMBER obser- 
vations. This favorable outcome of the analysis implies that the 
general morphology of the asymmetries in the disc of £ Tau are 
correctly reproduced by the model, at least in the regions of the 
disc where most of the Ha and Bry radiation originates. In ad- 
dition, the fit of the observed phase difference between the V/R 
cycles of Ha and Bry strongly supports the theoretical prediction 
that the density wave has a helical shape. 

However, two weak points of the model stand out from our 
analysis. Our results did not reproduce the phase of the V/R cy- 
cle of the Brl5 line, even though the amplitude was reproduced. 
Also, the model predicts large variations in polarization across 
one V/R cycle that are not observed. Both the Brl5 line and the 
polarization share a common property: both are formed in the 
very inner part of the envelope within a few stellar radii. 

Considering the above discussion, a pattern emerges. 
Observables that trace the larger scale of the disc (Ha and Bry) 
are reproduced by the model, whereas observables that map the 
inner disc (Brl5 and polarization) are not. This suggests that the 
global oscillation model fails to reproduce the shape and/or am- 
plitude of the density wave in the inner disc, but it does repro- 
duce the general features of the outer disc. 

We must bear in mind that there are at least two alternative 
explanations as to why the model does not fit the Brl5 V/R cy- 
cle. First, the site for line formation, as calculated by HDUST 
may be wrong. If the size of the Br 15 line emitting region were 
larger in our models, then the phase difference would have been 
smaller. Second, there are other physical processes that affect 
the structure of the inner disc, and these may also contribute to 
this discrepancy. For example, the non-isothermal effects on the 



Carciofi et al.: Testing the 2D Global Disc Oscillation Model 



13 



disc surface density calculated by ICarciofi & Biorkman (2008). 
Also, the fact that the disc is unlikely to be fed by a smooth 
mass transfer from the st ar, but instead by a series of outbursts 
(e.g. [Rivinius et al.L[l998l) . In this case the steady-state viscous 
decretion disc would not be a good representation for the inner 
disc. 

The simultaneous constraints of both the period and mode 
confinement, together with a sophisticated model, allow us to 
place more stringent limits on the V/R parameters than the pre- 
vious analyzes (e.g.. lVakili et al.L [l998). However, there are still 
uncertainties in some parameters. This is particularly true for 
the adopted hypothetical radiative force used to adjust the oscil- 
lation period. There is also uncertainty arising from the fact that 
we cannot constrain each parameter independently. For exam- 
ple, a larger value of rotation parameter k2(Q./£l c rit) 2 makes the 
mode more confined, whereas the mode becomes less confined 
for a larger value of viscosity parameter a. Thus, a similar good 
fit can be obtained for different sets of ^(^V^crit) 2 and a. In or- 
der to remove such a degeneracy, a more precise determination 
of the basic stellar parameters must be made. 

5. Conclusions 

We have successfully completed a very detailed modeling of the 
Be star £ Tau, using a very large set of observations that were 
reported in the first paper of this series, f Tau provides a unique 
opportunity for testing theoretical ideas about disc formation and 
structure and about the origin of V/R variations. First, the disc is 
known to have its global properties constant since at least 1992. 
Second, the presence of the binary companion makes it reason- 
able to assume that the disc is truncated at the tidal radius of 
the system. Finally, the inclination angle of the disc is very well 
constrained. 

We have assumed that the unperturbed disc structure is that 
of a steady-state (i.e., constant decretion rate) viscous Keplerian 
disc. We further assume that the disc is in vertical hydrostatic 
equilibrium, and to describe the density perturbations we have 
adopted the global osc illation model of lOkazakil (119911) and 
iPapaloizou et all (1 19921) . 

We employed the computer code HDUST of 
ICarciofi & Biorkmanl {2006) to solve the coupled prob- 
lems of radiative transfer, radiative equilibrium and statistical 
equilibrium and we fit, in a simultaneous and self-consistent 
way, all of the available observations, including the recent 
VLTI/ AMBER observations made in December, 2006. 

The main conclusions of our data analysis are 

- We have obtained an excellent fit of the global, time averaged 
properties of the disc. From the theoretical point of view, this 
represents strong evidence that the viscous decretion disc 
scenario is the mechanism responsible for the formation of 
the CS disc of f Tau. 

- The agreement between our predictions and the 
VLTI/AMBER observations and the acceptable fit of 
the V/R cycle of the Ha and Bry lines represent the 
first fairly comprehensive quantitative test of the global 
oscillation scenario. By fitting the spatially resolved inter- 
ferometric observations we have shown that our model can 
reproduce the general features of the outer disc quite well. 
In particular, this demonstrates that the density wave is a 
spiral, as predicted. 

- The model failed to reproduce the V/R cycle of the Br 15 
line and predicted a large variation in the polarization dur- 
ing a V/R cycle, which is not observed. Since both the Brl5 



line and the polarization comes from the inner disc within 
a few stellar radii, this suggests that the global oscillation 
model does not predict the correct geometry in this region. 
However, other possibilities for the discrepancy between the 
model and the observations exist, including non-isothermal 
effects in the inner disc structure and perturbations of the in- 
ner disc by outbursts. A detailed analysis necessary to test 
both ideas will be left to future work. 
- With the aid of the VLTI/AMBER observations we have, for 
the first time, completely determined the geometrical prop- 
erties of a Be star system. Non-spatially resolved determina- 
tions of the inclination angle (as, for instance, obtained from 
polarization studies) cannot distinguish between i < 90° 
(i.e., a disc whose northern face is towards the Earth) and 
i > 90° (i.e., a disc whose southern face is towards the Earth). 
Our detailed fit of the interferometric data allowed us to de- 
termine that the inclination angle of f Tau is 95°, i.e., the disc 
is seen almost edge-on and its southern face is towards us. 

Recently, there has been importan t progress in the theory of 
global disc oscillations. Ogilvie(2008) found that the one-armed 
oscillation modes are naturally confined in discs larger than sev- 
eral tens of stellar radii, when the 3D struct ure of the oscill ation 
mode is taken into account. According to Ogilvie (2008j), the 
vertical oscillation structure should not be neglected, because 
the vertical gravitational acceleration around an elliptical orbit 
of each gas particle in a disc perturbed by a one-armed oscilla- 
tion mode excites an oscillatory vertical motion. 

Given that R t < 20 R e , the one-armed modes are unlikely to 
be confined in the ( Tau disc even if this effect is taken into ac- 
count. The amplitude of the oscillation, however, is expected to 
be significantly influenced in the outer part of the disc where the 
major contribution to the Balmer emission arises . In the next pa- 
per of this series, we will apply the treatment by Ogilvie (2008) 
to the £ Tau disc, in order to study the effect of the vertical os- 
cillations on this star and see if the model parameters can be 
well constrained and if these parameters agree with the values 
obtained by this study. 
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